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Abstract 

We consider the lateral difFusion of a protein interacting with the curvature of the mem- 
brane. The interaction energy is minimized if the particle is at a membrane position with 
a certain curvature that agrees with the spontaneous curvature of the particle. We employ 
stochastic simulations that take into account both the thermal fluctuations of the membrane 
and the diffusive behavior of the particle. In this study we neglect the influence of the particle 
on the membrane dynamics, thus the membrane dynamics agrees with that of a freely fluc- 
tuating membrane. Overall, we find that this curvature-coupling substantially enhances the 
diffusion coefficient. We compare the ratio of the projected or measured diffusion coefficient 
and the free intramembrane diffusion coefficient, which is a parameter of the simulations, with 
analytical results that rely on several approximations. We find that the simulations always 
lead to a somewhat smaller diffusion coefficient than our analytical approach. A detailed 
study of the correlations of the forces acting on the particle indicates that the diffusing inclu- 
sion tries to follow favorable positions on the membrane, such that forces along the trajectory 
are on average smaller than they would be for random particle positions. 

1 Introduction 

During the last decade it has become more and more apparent that lateral difFusion of proteins in 
membranes plays a crucial role in cellular functioning.^'^ Therefore, a whole range of experimental 
techniques has been developed, which is constantly being improved in order to determine accurate 
values of lateral protein diffusion coefficients.'^''^ The most important methods include fluorescence 
recovery after photo bleaching, ^'^ fluorescence correlation spectroscopy,^'^ or single particle track- 
jj-^g 9,10 -^jjjig a large amount of data has been collected with these techniques the interpretation 
of results always depends on reliable models for the diffusive process. In some situations, like 
restricted difFusion due to corrals, ^''^^^ a certain, often rather crude, qualitative explanation is 
easily found, in other situations, however, this is by no means the case and a reliable quantitative 
interpretation can only be achieved, if corresponding theoretical calculations or simulations are 
performed. 

Only recently an increased interest in lateral diffusion has emerged from a theoretical viewpoint. 
In order to compare theoretical results with experiments it is necessary to take into account 
various aspects of the particular system. These include the nature of the membrane the particle 
is diffusing in, the properties of the diffusing particle, or the experimental method with which 
diffusion coefBcients are determined. A very important aspect in both theoretical calculations and 
the analysis of experimental results is that the membrane must not be regarded as a flat plane but 
is often structured such that regions with higher and lower curvatures appear. For example, this 
must be accounted for in the study of diffusion in the membranes of the endoplasmic reticulum. 
Neglecting the influence of the membrane shape leads to considerable errors in the determination 
of diffusion coefficients.^'^ Several analytical and simulational studies have been performed that 
regard diffusion on various fixed curved surfaces. ^"^^^^ 

But even if a membrane appears to be flat on average, it is subject to thermal fluctuations 
that lead to rapid shape changes around the flat configuration. ^^'^^ Neglecting the influence 
of the membrane on the movement of the particle these fluctuations that depend on properties 
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Figure 1: Sketch of the curvature-coupled model for an inclusion at position R, with area TrOp, 
a bending rigidity m and a spontaneous curvature Cp in a model membrane with the bending 
rigidity k and the effective surface tension a. The membrane is mapped to a two-dimensional 
N X N lattice with the lattice spacing £ of which only a one-dimensional cut is shown here. The 
lateral dimension of the system is L = iV x ^. In the simulations, the area of the inclusion is 
determined by the weighting fimction G(r — R), which is set to be a Gaussian. 

like bending rigidity, surface tension, proximity to substrates or other membranes, etc., have an 
influence on the measured values for lateral diffusion coefficients because experiments usually 
regard the path of the inclusion projected on a flat reference plane instead of the actual path 
along the membrane. Compared to the intramembranc diffusion coefficient the measured diffusion 
coefficient will be the smaller the stronger the fluctuations. This was initially pointed out by 
Gustafsson and Hallc;^'' the quantitative evaluation of this effect for free intramembranc diffusion 
was performed more recently in independent work by us and two other groups using analytical 
calculations^^'^^ and simulations. ^"'^^ 

All studies mentioned in the last two paragraphs take into account the influence of the shape 
of the membrane on measured diffusion coefficients, but otherwise neglect any interaction between 
membrane and protein. Considering that a protein also has certain physical properties it must 
be assumed that interactions between the protein and the membrane exist that influence the 
diffusion coefficient. This assumption is corroborated by experimental findings: for example after 
photoactivation of bactcriorhodopsin (BR) in model membranes the lateral diffusion coefficient 
is reduced by a factor of five.^^ This reduction is attributed to oligomerization of BR upon 
activation caused by structural changes that influence protein- lipid interactions. There are many 
other experimental examples that indicate that the interactions between membrane and protein 
have a signiflcant influence on lateral diffusion, see for example refs.^*'-^^ An important property 
of a membrane compared to a flat surface is the curvature. A variety of studies mainly using 
particle based simulations are concerned with the influence of inclusions with a certain intrinsic 
curvature on membrane shape and lateral diffusion. "^^"^^ In earlier work we calculated effective 
diffusion coefficients for particles with a bending rigidity and a spontaneous curvature. These 
calculations revealed that the additional interaction that tries to move the particle to positions on 
the membrane where the curvature agrees with the particle's spontaneous curvature, leads to an 
increase in the diffusion cocfRcient. 

Our previous work on curvature-coupled diffusion effectively describes diffusion of a point- 
like particle and relies on several approximations. ^'^'^^ One of these, the so-called pre-averaging 
approximation, assumes that membrane fluctuations on all possible length scales have a much 
shorter relaxation time than the time it takes the particle to diffuse the corresponding distance. 
In this paper we introduce a scheme to simulate the diffusion of a particle with a certain extension, 
a bending rigidity and a spontaneous curvature in a thermally fluctuating membrane. A sketch of 
the considered system with the relevant physical parameters is given in fig. [TJ 

The present method is no longer restricted to certain relative timescales of diffusion and mem- 
brane fluctuations. The additional energy of the inclusion is introduced by replacing a patch of 
membrane that is described by the Helfrich Hamiltonian with a new Helfrich-likc term with a 
different bending rigidity and a spontaneous curvature. The extension of the particle is modeled 
by a Gaussian weighting function in order to have a smooth crossover from the bare membrane to 
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the particle. Obviously the system is dominated by two dynamic processes: the shape fluctuations 
of the membrane and the particle diffusion. Two Langevin-equations, one for the membrane, the 
other for the particle, are derived from the energy of the system. Our simulation scheme consists of 
the numerical integration in time of these coupled equations. Apart from performing simulations 
we also analytically evaluate the coupled dynamic equations by use of a perturbation theory that 
neglects the influence of the particle on the membrane and assumes that membrane relaxation 
times are much smaller than corresponding diffusive time scales. In order to compare with these 
analytical calculations and to reduce the computational effort, we restrict our simulations in the 
current study such that the membrane dynamics is also not influenced by the diffusing particle. 
The influence on membrane movement will be considered in future work. The main quantity of 
interest is the ratio of the curvature-coupled and the intramembrane diffusion coefficient as a func- 
tion of the membrane parameters bending rigidity and surface tension. The latter coefficient is a 
parameter of our scheme and resembles the free diffusion coefficient of the particle if no additional 
force were acting on it. The application of both our approaches shows that curvature-coupling 
leads to increased diffusion. However, the comparison of our analytic calculations with the sim- 
ulation results reveals a systematic difference. In order to gain insight into the reason for these 
discrepancies we study force correlation functions that are the main contributions to the diffusion 
constant. 

The paper is organized as follows: In the following section we explain the model for the 
membrane dynamics and the Langevin-cquation for the inclusion. In sec. [31 the method and 
the approximations of the analytical calculation of the curvature-coupled diffusion coefficient are 
presented while sec. |4] introduces the used simulation scheme. Sec. [5] discusses the choice of 
parameters used in both the calculations and the simulations. The presentation of the results 
in sec. |6]is followed by a detailed discussion in sec. [71 why simulations lead to smaller diffusion 
constants and a possible interpretation of our findings. The paper finishes with some conclusions 
and an outlook for future work. 



2.1 Membrane dynamics 

We consider a model membrane in a fluid environment. The membrane is given in Monge- 
reprcscntation where r = (a;, y)^ is the position in the (x, ?;)-planc with the deviation /i(r, t) 
out of this plane. For such a membrane Helfrich'^^ derived the free energy that to lowest order has 
the following form in the Monge-gauge'*'' 



with the bending rigidity k, the effective surface tension a and the area in the [x, j/)-plane. 
The dynamics of a membrane is given by^^ 



with the Onsager coefficient A(r — r') that takes into account hydrodynamic interactions with the 
fluid background. Using the Fourier-transformation 
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one obtains ([2]) in Fourier-space 
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with E{k) = Kk^ + crfc^ and A(fc) = 1/ {Arjk) which is the Fourier-transformed Onsager coefficient 
for a free membrane in a fluid with viscosity rj.^^ The stochastic force <^(k, t) obeys the fluctuation- 
dissipation theorem 



(e(k,t)) = 0, 

(c(k,or(k',t')) 



2A{k)—6{t-t')S^,^, 



(6) 
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where /? = {ki,T) ^ is the inverse temperature. Later on, we need the correlations of h(k, t) which 
are given by 
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2.2 Diffusion 

Now we place an inclusion into the membrane that diffuses freely along the membrane. The 
dynamics of the inclusion may be described by a Fokker-Planck-equation (FP-eq.). However, since 
the diffusive motion takes place on a curved surface the Laplace-operator needs to be replaced by 
the Laplacc-Beltrami-opcrator. This leads to a new FP-eq.^"''^^ 



dtP{r, t)=DY, d^^g^'d—P (r, t) , 



(10) 



with the diffusion coefficient D, the metric g and the inverse metric tensor g^^ . In the Monge 
gauge the metric has the form g = 1 + h'^ + hy, while the inverse metric tensor is 
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(11) 



The subscripts denote partial derivatives, e.g. = dh/dx. The probability P{r,t) of finding the 
projection of the inclusion at a position r is normalized to JdrP{r,t) — 1. In the simulations 
we make use of a Langcvin-equation to describe the motion of the projected particle position 
R{t) = {X{t),Y{t))'^ . Using the above FP-eq. a projected Langevin-eq. is derived within the 
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Stratonovich calculus:^^'"'^ 
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The upper case subscripts express that the partial derivatives at the particle position R(t) have 
to be used. The stochastic force ^ has zero mean and is delta-correlated: 



(G> = 0, 



(13) 
(14) 



Equation (|12p comprises a drift that is caused by the membrane curvature and diffusive terms. 
The consequences of such a drift term for a freely diffusing inclusion have been introduced in ref.^^ 
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2.3 Curvature-coupled model 

The equations derived so far apply to a freely diffusing point-like inclusion. If one is interested in 
the diffusion of a more realistic inclusion, one has to take the physical parameters of the inclusion 
into account. First the inclusion has a non vanishing area which is set to Tra^. Furthermore, the 
inclusion has possibly its own bending rigidity m and maybe a spontaneous curvature Cp. As 
indicated by "into the membrane" the inclusion completely replaces the membrane at its position. 
To consider this in the free energy of the system, one has to add a new energy term for the 
inclusion and remove the part of the membrane which is replaced. The additional term caused by 
the inclusion leads to the new free energy 

n = no + ni, (15) 

where 



Hi [h{r,t),R{t)] = J drG(r-R)x 



^(V?Mr,i)~C,)^-^(V?Mr,t))'] (16) 



is the correction to Helfrich's free energy Hq. G'(r — R) is a weighting function for the extension 
of the particle that we set to be a Gaussian such that the crossover from particle to membrane is 
smooth. Taking into account the area constraint JdrG(r — R) = Tra^ it is given by 

G(r-R)=exp|-^^— (17) 

The altered free energy ()16p induces additional forces on the inclusion and the membrane. The 
membrane dynamics is obtained by replacing TCq in ([2|) with Ti. such that 

SMr, t) ^ - /*'A(r-rO (^ + ^) + (18) 

The forces that influence the diffusive behavior of the inclusion can be calculated by f = — VrTYi. 
Taking into account the curvature of the membrane in the force term, which needs to be added to 
the right hand side of (fT2|l . we get the complete equation of motion for the inclusion 

^tR^ (t) = ^tRproJ,^ ' f^'Y. 9^' ^.Hl, (19) 

j 

with the mobility fi that is related to the intramembrane diffusion coefficient D via the Einstein 
relation D = kbTfi. 

With eqs. p8|) and (flQl) the dynamics of the system is fully determined. Note, that these 
equations are coupled since the particle diffusion depends on the shape of the membrane via the 
partial derivatives of ft.(r, t) at the particle position, and the membrane dynamics on the position 
of the particle through the additional energy. 

3 Analytical calculations 

In order to calculate a new curvature-coupling affected diffusion coefficient defined as 

^cc ^lim \ ' 20 

one has to determine the mean square displacement 

t t 

(AR2(i)) ^JdrJ dr' (9,R(t) • dr'Rir')) , (21) 
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by integrating eq. ([TO]) in time and performing the thermal average. Since the expHcit calculation 
of the mean square displacement using the exact equation of motion (jl9p and the full membrane 
dynamics (|18p is not possible analytically, it is necessary to introduce several approximations. 
In order to simplify eq. (jl9p we perform a pre-averaging approximation. This approximation 



is applicable if for all modes the membrane relaxation times {A{k)E{k))~^ , see eq. (|9]), are con- 
siderably shorter than the time tt"^ /{Dk"^) it takes a particle to diffuse the distance given by the 
corresponding wave length of the mode. For typical experimental values for bending rigidity k, 
tension a, diffusion coefficients D, and system sizes L, this condition is very often fulfilled. If 
membrane fluctuations are "faster" than the diffusion of the particle it is assumed that the par- 
ticle only feels average membrane fluctuations. The applicability of this approximation for free 
lateral diffusion is discussed in ref.^^ In our current work the pre-averaging approximation results 
in the replacement of fi-g^-' in eq. (|19p by fiproj^ij where the projected free mobility is defined by 

We, furthermore, assume that the additional energy Tii caused by the insertion of a single 
particle is small, in order to justify a perturbation expansion to first order in the particle energy. 
The consequence of this approximation is that the dynamics of the membrane is not influenced by 
the presence of the inclusion. Thus the membrane dynamics is expressed by eq. ([5]). 

Another approximation needs to be employed so as to make analytical calculations possible. 
Inserting eq. (|19p into eq. ()2ip we see that the mean square displacement becomes a function 
of the height correlations {h{'R{t),t)h(R{t'),t')) . These correlations decay with increasing time 
difference |i — t'| due to two reasons: during the time interval the membrane shape changes and 
the particle position advances. Since we assume diffusion to be much slower than membrane shape 
changes we neglect the effect caused by the particle movement. 

Note, that this approximative analytical calculation cannot include a possible correlation be- 
tween the particle position and the membrane shape in the vicinity of the inclusion. In other 
words, we assume a constant probability for finding the particle at any point in the system rela- 
tive to a given membrane configuration. This aspect becomes important when we compare with 
simulation results, as will be discussed in sec. [T) 

Using the inverse Fourier-transform given in eq. Q and applying the previously explained 
approximations we find for the mean square displacement 

(AR2 (<)) = ADprojt+ 
t t 

" " k k' 



(22) 
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with the hight correlation function given in cq. Q and Dp 



Inserting the resulting equation for (AR^ (t)) into (j20p and performing the long time limit one 
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(23) 



In this equation it is interesting that there is no need to set a cut-off for the wavenumber k since 
the exponential function cxpj — fc2ap/2} damps higher k values. In ref.^^ a similar calculation 
for the curvature coupled diffusion coefficient is performed. There the area function has the form 
TTCpJ (r — R) and a cut-off is necessary. The resulting diffusion coefficient agrees with the above 
diffusion coefficient of equation (|23p in the limit of vanishing variance of the Gaussian. 
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4 Simulation method 



To probe the applicability of our analytical calculations that depend on several approximations we 
set up simulations that numerically integrate the coupled equations of motion for the membrane 
and the diffusing particle. We use a square, periodic lattice with N x N lattice points and 
the lattice spacing £ to map a model membrane with size L = N x £, sec fig [TJ To simulate 
the shape fluctuations of the membrane wc numerically integrate the appropriate equation of 
motion. In order to reduce the computational effort and to compare with the analytical calculations 
introduced in the previous section we evaluate the unperturbed Langevin-equation ([5|) discretely 
in time. These calculations are performed in Fourier-space since the equations of motion for the 
height function modes /i(k, t) decouple. Due to the periodic boundary conditions the wave vectors 
are of the form k = 2'K{l,n)/L with I and n being integers. Since the hight function /i(r) is a 
real function defined on a iV x iV lattice the relation ft,(k, t) = k, applies leading to the 
restriction -N/2 <l,n^ N/2. 

Regarding the fluctuation-dissipation theorem ([7]) it is obvious that fluctuations of h{k = 0,t) 
would diverge for the Onsagcr-cocfficient A(k ^ 0) ^ oo. The dynamics of the height function 
mode ft,(k = 0,t) corresponds to the center of mass movement of the whole membrane. Due to 
the irrelevance of this movement in the determination of the lateral diffusion coefficient wc set 
A(k = 0) = and keep /i(k = 0, = fixed at all times. 

In order to choose an appropriate discrete time step At that ensures that numerical errors are 
small it is necessary to point out that the largest k-vcctor possible kmax = V2-kN/L in the given 
lattice determines the smallest time scale of the membrane as can be seen in eq. The used 
time step At in the simulation should be smaller than this smallest time scale. 

The dynamics of the inclusion is given by a discrete version of cq. (fT5|) that is also numerically 
integrated in time. Since it is coupled to the membrane equation ^ via the derivatives hx, hy, 
etc., of the membrane configuration at the position of the inclusion, the temporal evolution during 
a discrete time step At consists of an update of the membrane shape and the particle position. In 
contrast to the membrane dynamics the motion of the inclusion is calculated in real space. The 
required derivatives of /i(r, t) are, therefore, determined in Fourier-space and then transformed to 
real space using routines of the FFTW- libraries. ^"^ We allow off-lattice diffusion for the inclusion 
which is necessary if the average time it takes a particle to diffuse a distance £ is much larger than 
the corresponding membrane relaxation time. Thus the derivatives at the position of the inclusion 
are determined by a distance weighted linear extrapolation of the four nearest lattice sites. In 
order to minimize the numerical error caused by not following the membrane surface correctly the 
displacement per time step is to be set to a small fraction of the lattice spacing £. With the above 
explained restriction for the time step in order to describe the membrane shape evolution with 
sufficient accuracy there are overall two conditions that must be fulfilled in the choice of At. 

Apart from obvious computational limits the determination of how long simulation runs should 
at least be is again dictated by two time scales. On the one hand the length of the simulation 
should be several times the longest relaxation time of the membrane, which is given by the smallest 
possible wave vector, to ensure that the membrane shape has passed through an adequate amount 
of independent configurations. On the other hand it is preferable that the inclusion has enough 
time on average to cover a distance of several lattice sites. 

A more detailed description of the simulation method is given in ref.,^^ where the corresponding 
scheme for free particle diffusion has been introduced. 

5 Parameters 

Before we present our results we will introduce the used parameters. These are the same for 
the analytical calculations and the simulations. As already mentioned in the description of the 
simulation method we use a discrete membrane with a lattice spacing of £ that we set to £ = lOnm. 
This choice reflects a compromise between the wish to simulate reasonably sized systems and the 
computational limitation in the number of lattice sites. A decrease in the lattice spacing for a 
constant system size, i.e. an increase in the number of lattice sites, introduces additional large k- 
values that contribute only weakly to membrane fluctuations, see cq. (j9l), or the curvature-coupled 
diffusion coefficient £ is one of the basic units. The others are the time given in seconds s 
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and the thermal energy that is (3^^ ~ 4.14 x lO^^^J at room temperature. All parameters of the 
system are given in units of s, and (3. For the determination of the membrane parameters we 
look at typical experiments and extract a range of 5 to 50 for the bending rigidity k/3 and 10~^ 
to 10~^J/m^ for the effective tension a that corresponds to 10~^ to 10^^ for a(3£'^. At around 
af3£'^ = 10 rupture of the membrane occurs. In the experiments it is, furthermore, common to use 
water as a surrounding medium with a viscosity of 77 = 10~'^kg/ms = 2.47 x 10~^s/ [(3£'^). The last 
parameter of the membrane is the size L that is related to the number of lattice points N in each 
direction via L = N x £. Since a sufficient number of wave vectors k are considered for a 50 x 50 
lattice we use a system size of L = 50£ that corresponds to 0.5/im. For the parameters of the 
inclusion we choose Cp£ = 2, mf3 = 2k(3 and Up — 1£ as in our previous calculations.^^ The bare 
diffusion coefficient D, however, has to be chosen carefully. Since the analytical calculations rely 
on a pre-averaging approximation we have to set D sufficiently small. Therefore, the time scale of 
the diffusion td — tt^/ {Dk^^^^ must be much longer than the largest membrane time. This time 
is given by Tmemb,max = / {^f^min + ^^min) with the absolutc value of the smallest wave vector 
kmin = 211 / L which corresponds to the longest length L in the system. The comparison of the 
two time scales leads to L» < tt^ {nkmm + <^Kiin) I (4*7) ~ W'£^ js. We choose = 5 x 10''^^/s. 
This choice is good for all values of k/3 and in the selected range. For the total length of the 
simulations one has to keep in mind that it has to be several times longer than the time scale of 
the longest wave vector. The time step At however has to be so small that it is smaller than the 
shortest membrane time T,„em(,_m.i„ and that the inclusion diffuses only a short distance. We set 
the total length to 1ms and At = 10~^s such that each simulation run comprises 10® time steps. 
This choice of the time step At is applicable for small bending rigidities k/3. 

Note, that an increase of k/3 and cr/3£^ leads to smaller time scales Tmemb of the membrane. 
In order not to increase computing time we keep At fixed for all considered k/3 and a/3£'^. This, 
however, will lead to a slight increase in numerical errors for the modes ft,(k, t) with very large 
wave vectors in membranes with large Kf3 and a/3£'^. Since fluctuations of these modes, see eq. ([9|). 
are rather small these errors are negligible in the determination of Dec- 



6 Results 

The analytically derived curvature-coupled diffusion coefficient D^c is calculated by numerical 
summation of equation (|23p using Mathematica for the whole range of parameters given in the 
previous section. The ratio Dcc/D of the curvature-coupled Dec and the free intramembrane 
diffusion coefficient D is plotted as a function of the bending rigidity k/3 and the effective surface 
tension (t/3£^ in fig. [2] The ratio increases for brighter colors. In the figure one can see that for small 
al3£'^ the ratio increases very strongly for k(3 < 5. With a further increase of the bending rigidity 
k/3 the ratio reaches a plateau. Here the strengthening of the forces f caused by an increasing 
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Figure 3: Mean square displacement determined from simulation data as a function of time for 
Kf3 — 5 and cr/3£^ — 0. Plotted with the data points are the two fits for the diffusion coefficient at 
the beginning and the end of the simulations. Furthermore, the expected mean square displacement 
for bare free diffusion is displayed. 




Figure 4: Ratio Dec/ D of the curvature-coupled to the intramembrane diffusion coefficient as a 
function of bending rigidity Kf3 for fixed surface tension a(3i'^ = 1 x 10^^. Most of the simulated 
data points are smaller than the analytical curve determined by eq. (j23p . 



bending rigidity of the inclusion to/3 is compensated by the fact that thermal fluctuations become 
weaker for increasing The increase of (t/3£^ by about three orders of magnitude leads to no 
significant effect. However, for even larger values of cr/3£^ the ratio decreases fast to one for small 
K/3. This happens since a large surface tension af3i'^ also damps the thermal fluctuations of the 
membrane and, therefore, the additional force on the inclusion is small. Increasing k(3 the ratio 
also increases for large tr/?^^ but does not reach the same hight as for small tensions. In this case 
the increase of Dec/ D for high values of K.f3 is also compensated by the damping caused by the 
surface tension. 

Overall, our calculations show that the curvature-coupling, which leads to an additional force 
f on the inclusion, enhances the inclusion's diffusion rate in the investigated parameter range. It 
is noteworthy that the ratio is always bigger than one despite the fact that the projection alone 
would lead to a ratio smaller than one.^^ 

Results from variations of the other parameters, like L, Cp, etc., are not plotted but the effect 
on Dec can easily be obtained from eq. ([23|) . However, one has to keep in mind that the effect of the 
inclusion on the membrane has to be small in order for the perturbation theory to be applicable. 

In the simulations, we also use the unperturbed membrane equation of motion ([5]) and the same 
parameter sets as for the numerical summation just discussed. As we are interested in results in 
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Figure 5: Ratio Dcc/D as a function of cr/3£^ for the fixed bending rigidity n(3 = 10. The simulated 
values (+) follow the analytical curve (line) of eq. ((23|) only qualitatively. 

thermal equilibrium each membrane starts in a random configuration and has 1 ms to equilibrate 
before the inclusion is placed in its center. 1 ms is about five times the longest membrane relaxation 
time so we can be sure that the membrane is in thermal equilibrium. Since we obtain only one 
particle trajectory per independent simulation run we average over 500 simulations with the same 
set of parameters to get the mean square displacement <|AR^(t)). An example for ^AR^(t)) as 
a function of time is plotted in figure [3l The slope of the resulting straight line at late times 
corresponds to Dec, see eq. (|20|) . The determined values for Dcc/D from the simulations (+) are 
plotted with the analytical curve (line) in fig. [4] as a function of k(3 for a constant a^l"^ and in 
fig. E] as a function of apf' for a constant bending rigidity k/3. We see that the simulations follow 
qualitatively the analytical curve but the values are about 10% smaller than expected. 

7 Discussion 

Both of our approaches demonstrate that curvature- coupling enhances diffusion. This result is 
plausible for the following reason. Due to the membrane fluctuations the positions that are favor- 
able for the diffusing particle are constantly changing. Thus the particle is subject to changing 
forces leading to enhanced movement of the particle, which in turn leads to a higher diffusion co- 
efficient. The resulting enhanced diffusion coefficient is caused by forces, which are still thermal. 
Therefore, the system is still in equilibrium and the ffuctuation-dissipation-theorem is applica- 
ble, such that an increased effective mobility or a reduced effective friction of the particle can be 
determined. 

On the quantitative side, the analysis of the simulation data reveals a diffusion coefficient that 
is about 10% smaller than we expect from the analytical calculations. In these calculations several 
approximations that we have explained in sec. [3] are applied to calculate the mean square displace- 
ment (AR2(t)) dH])- The dominant contribution to Dec is (f(r, R(r)) • f (t -h Ar, R(r At))) 
with the force f (r, R(r)) = — VRHi[/i(r,T),R(r)], see eqs. (PT|) . Hence we investigate this 

force correlation function. 

First we consider the averaged quadratic force for At = along the trajectory R(t) of the 
inclusion to see whether differences in the strength of the correlations occur. Then we will regard 
the dependence of the correlations on the time interval At. To examine the strength of the 
correlations for At = the quadratic force acting on the inclusion is determined at each time 
step of the simulation and then averaged over all values. For several sets of parameters we receive 
the values that are plotted in fig. [6] as a function of cr/3£^ for a constant k/3 and in fig. [7] as a 
function of k/3 for a constant a(3^'^ . The values resulting from the analytical calculations are also 
plotted in these figures. We observe a difference between analytical and simulation results that 
is on the same order of magnitude as the difference in the diffusion coefficients. In the analytical 
calculations we assume that the probability of finding the inclusion at a particular position is the 
same for any point of the membrane. If we calculate, using the simulation data, the mean square 
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Figure 6: Average over the squared force on the inchision for equal times (f^ (r, R(t))) as a 
function of the effective tension a(3i'^ for the fixed bending rigidity k/3 ~ 5. The figure shows a 
good agreement between the analytical points (<)) and those achieved by averaging over a fixed 
point on the lattice (x). The averaging along the trajectory (+) leads to smaller values. 
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Figure 7: Average over the squared force on the inclusion for equal times (f^(r, R(r))) as a 
function of the bending rigidity k(3 for the fixed surface tension a^i"^ = 0. The figure shows a 
good agreement between the analytical points (<)) and those achieved by averaging over a fixed 
point on the lattice (x). The averaging along the trajectory (+) leads to smaller values. 
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of the force the inclusion would be exposed to if it were fixed at some arbitrary position on the 
membrane we obtain values very close to the mean squared force resulting from the analytical 
calculations. These values are also plotted in figs. [S] and [71 The differences in fig. [B] between the 
analytical calculations and the average for a fixed point are caused by the time step At of the 
simulations that induces bigger numerical errors for higher values of k/3 as previously explained in 
section [5] Overall, averaging over the whole lattice in the calculations leads to seemingly higher 
forces than along the actual particle trajectory. A possible explanation for this reduced force is 
that most of the time the inclusion is close to a local minimum of the free energy. The extrema 
of the energy are created by the membrane shape, which is constantly changing due to thermal 
fluctuations. Therefore, the positions of the extrema will also move along the membrane. As the 
negative gradient of the energy is always pointing towards the nearest local minimum the inclusion 
will predominantly move in the direction of the nearest local minimum. For a fast enough diffusion 
rate the inclusion is capable of following a local minimum. 

Due to the fact that for each simulation run a new thermally equilibrated membrane is used and 
the particle is always placed in the center of the membrane, it is very unlikely that the inclusion 
is initially close to an energy minimum. Therefore, the inclusion is exposed to higher forces at the 
beginning than at later times. Since higher forces go along with a higher diffusion rate we expect 
to observe two diffusion coefficients from the simulation data: a smaller one for late and a larger 
one for early times. Indeed such a behavior occurs as one can see in fig. [3] where the linear fits to 
the mean square displacement for early and late times arc plotted. In this example the crossover 
is at about 0.2ms. Comparing the resulting diffusion coefficients with the analytical values, the 
one determined at the beginning of the simulation agrees reasonably well with the analytically 
determined diffusion coefficient. The inclusion starts with a higher diffusion rate and then, after 
a variable time period, finds a local minimum, which it tries to follow. 

Now that we have found that the force acting on the inclusion is reduced in the simulations we 
consider the time correlations of the force (f (r, R(t)) • f (t + At, R(r + At))) . We expect to see a 
reduction of the correlations caused by the inclusion following an energy minimum but in addition 
the time dependency will be influenced by the movement of the particle position R(t). The 
obtained correlation function along the trajectory is plotted in fig. [8] together with the analytical 
result for nf] = 5 and vanishing tension. For a better representation the functions are normalised 
to one and plotted for small Ar in the inset of the figure. It is obvious that in addition to the 
altered start values the time dependence is also different. The decay of the correlations along the 
trajectory is slower than for the analytical curve. To demonstrate that this altered decay is caused 
by the motion of the inclusion we choose five fixed points of the lattice and determine, from the 
simulation data, the time correlation function of the force that would act on the inclusion if it were 
fixed at these points. The average over these points is also plotted in fig. [5] and agrees well with the 
analytical curve. The altered decay along the trajectory is an effect caused by the motion of the 
particle and indicates that the force correlations are stronger along the trajectory. This fact does 
not only lead to a slightly higher diffusion rate but also shows that the forces for a series of time 
steps point in similar directions which corroborates our interpretation that the inclusion tries to 
follow a local minimum. Although a larger decay time of correlations leads to enhanced diffusion 
the observed diffusion rate is smaller than the analytically calculated value since the effect that 
forces close to energy minima arc reduced dominates. 

8 Conclusions 

In this paper, we have derived a model for the interaction of an inclusion with a model membrane 
and have investigated the influence of this interaction on the diffusion of the inclusion. In the model 
the inclusion is a physical object with an area, a spontaneous curvature and a bending rigidity. 
In analogy to Helfrich's free energy a new free energy is derived and with this coupled equations 
of motion for the inclusion and the membrane dynamics. Using these stochastic equations and 
employing several approximations we calculate the curvature-coupled diffusion coefficient Dec- To 
assess the quality of the approximations in this analytical calculation that assumes a weak pertur- 
bation of the membrane we set up simulations. These simulations that numerically integrate the 
coupled equations of motion for the membrane and the inclusion are also based on the unperturbed 
membrane equation. 
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Figure 8: Time correlation function of the additional force on the inclusion 
(f (r, R(t)) • f (t + Ar, R(t + At))) for the bending rigidity k/3 = 5 and the effective ten- 
sion apP = as a function of the time difference Ar. In the inset the first 10~^s of the time 
correlation function are plotted scaled to one for At = 0. A good agreement is found for the 
analytical curve (line) and the average over five fixed lattice points (x). The time dependence of 
the force correlations along the trajectory (+) is different from the analytical result. 



Both, the simulations and our analytical approach, clearly display that the additional force 
on the inclusion caused by the interaction between particle and membrane leads to a significant 
increase in the diffusion coefficient compared to bare intramembrane diffusion. However, com- 
paring the analytical calculations with the simulation results one finds that the curvature-coupled 
diffusion coefficient in the simulations is about 10% smaller than analytically expected, but follows 
the behavior qualitatively. A closer look at the forces on the inclusion shows that the averaged 
forces along the trajectory of the inclusion are smaller than the averaged notional forces acting at 
any arbitrary fixed point of the lattice. The average over the latter forces, however, has a good 
agreement with the values from the analytical calculations. Since smaller forces correspond to 
local extrema of the free energy, and the gradient, i.e. the force, always points to the nearest local 
minimum it is likely that the inclusion is in the vicinity of such a local minimum of the energy 
most of the time. If the mobility of the inclusion is big enough the inclusion is able to follow such 
a local minimum. Another point for this interpretation is that wc place the inclusion in a new 
thermally equilibrated membrane for each simulation run. Hence, the inclusion is not necessarily 
close to a local minimum at the beginning and the diffusion rate should be faster for early times of 
the simulations. Considering the simulation data we find indeed that the diffusion for early times 
is about the same as the diffusion rate expected from the analytical calculation. After a short 
time the diffusion rate decreases and then remains constant. This corroborates the assumption 
that the inclusion meets a local minimum after some time and then tries to follow it. 

Our study leads to the conclusion that the analytical calculations provide qualitative values for 
the curvature-coupled diffusion coefficient for a given set of parameters. For a quantitative value 
of the curvature-coupled diffusion coefficient simulations are necessary that take the effects of the 
movement of the inclusion into account. 

In this paper we use the unperturbed membrane equation of motion. In ongoing work we 
are investigating the influence of the inclusion on the membrane dynamics and the diffusion by 
use of simulations that incorporate the perturbed membrane equation of motion. Then, the 
particle does not only adapt to the membrane but the membrane also adjusts to the particle. This 
will possibly lead to a further reduction of the diffusion coefficient compared to the analytical 
calculation. By comparing these results with those of the present paper it will be possible to 
determine the parameter range in which the perturbation may be neglected. We intend to also 
study the diffusion of several inclusions in a membrane in order to investigate possible cluster 
formation induced by membrane mediated interactions between the inclusions. Aside from this 
one may also be interested in other forms of interactions between the membrane and the inclusion 
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or additional interactions between the inclusions. 

The investigation of several possible interactions and the resulting effects on the diffusion co- 
efficient of inclusions as a function of membrane parameters may help to understand experimental 
data better and should finally lead to a deeper insight of the diffusion processes in biological 
membranes. 



14 



References 

1. Lippincott-Schwartz, J.; Snapp, E.; Kenworthy, A. Nat. Rev. Mol. Cell Bio. 2001, 2, 444-456. 

2. Marguct, D.; Lcnnc, P. F.; Rigncault, H.; He, H.-T. EMBO J. 2006, 25, 3446-3457. 

3. Lommcrsc, P. H. M.; Spaink, H. P.; Schmidt, T. BBA-Biomembranes 2004, 1664, 119-131. 

4. Chen, Y.; Lagerholm, B. C; Yang, B.; Jacobson, K. Methods 2006, 39, 147-153. 

5. Reits, E. A. J.; Neefjes, J. J. Nat. Cell Biol. 2001, 3, E145-E147. 

6. Lippincott-Schwartz, J.; Altan-Bonnet, N.; Patterson, G. H. Nat. Cell Biol. 2003, 5, S7-S14. 

7. Kohl, T.; SchwiUe, P. Adv. Biochem. Eng. Biot. 2005, 95, 107-142. 

8. Thompson, N. L.; Lieto, A. M.; AUen, N. W. Curr. Opin. Struct. Biol. 2002, 12, 634-641. 

9. Saxton, M. Annu. Rev. Biophys. Biomol. Struct. 1997, 26, 373-399. 

10. Kusumi, A.; Nakada, C; Ritchie, K.; Murasc, K.; Suzuki, K.; Murakoshi, H.; Kasai, R. S.; 
Kondo, J.; Fujiwara, T. Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 351-378. 

11. Tomishige, M.; Sako, Y.; Kusumi, A. J. Cell Biol. 1998, 142, 989-1000. 

12. Edidin, M.; Kuo, S. C; Sheetz, M. P. Science 1991, 254, 1379-1382. 

13. Sbalzarini, I. F.; Hayer, A.; Hclcnius, A.; Koumoutsakos, P. Biophys. J. 2006, 90, 878-885. 

14. Aizenbud, B. M.; Gcrshon, N. D. Biophys. J. 1982, 38, 287-293. 

15. Holyst, R.; Plewczyhski, D.; Aksimcntiev, A. Phys. Rev. E 1999, 60, 302-307. 

16. Plewczyhski, D.; Holyst, R. J. Chem. Phys. 2000, 113, 9920-9929. 

17. Faraudo, J. J. Chem. Phys. 2002, 116, 5831-5841. 

18. King, M. R. J. Theor. Biol. 2007, 227, 323-326. 

19. Yoshigaki, T. Phys. Rev. E 2007, 75, 041901. 

20. Naji, A.; Brown, F. L. H. J. Chem. Phys. 2007, 126, 235103. 

21. Brochard, F.; Lennon, J. F. J. Phys. (Pans) 1975, 11, 1035-1047. 

22. Seifert, U. Adv. in Phys. 1997, 46, 13-137. 

23. Gustafsson, S.; Halle, B. J. Chem. Phys. 1997, 106, 1880-1887. 

24. Reister, E.; Seifert, U. Europhys. Lett. 2005, 71, 859-865. 

25. Gov, N. S. Phys. Rev. E 2006, 73, 041918. 

26. Reister-Gottfried, E.; Leitenberger, S.; Seifert, U. Phys. Rev. E 20 07, 75, 011908. 

27. Kahya, N.; Wiersma, D.; Poolman, B.; Hoekstra, D. J. 5io/. C/iem. 2002, 277, 39304-39311. 

28. Vereb, G.; SzoUosi, J.; Matko, J.; Nagy, P.; Farkas, T.; Vfgh, L.; Matyus, L.; Wald- 
mann, T. A.; Damjanovich, S. Proc. Natl. Acad. Sci. USA 2003, 100, 8053-8058. 

29. Forstner, M.; Yee, G.; Parikh, A.; Groves, J. J. Am. Chem. Soc. 2006, 128, 15221-15227. 

30. Weikl, T. R.; Kozlov, M. M.; Helfrich, W. Phys. Rev. E 1998, 57, 6988-6995. 

31. Weikl, T. R. Europhys. Lett. 2001, 54, 547-553. 

32. Dommersnes, P. G.; Fournier, J.-B. Eur. Phys. J. B 1999, 12, 9-12. 



15 



33. Dommersnes, P. G.; Fournier, J.-B. Eur. Phys. J. E 2003, 11, 141-146. 

34. Fosnaric, M.; Iglic, A.; May, S. Phys. Rev. E 2006, 74, 051503. 

35. Cooke, I. R.; Dcscrno, M. Biophys. J. 2006, 91, 487-495. 

36. Reynwar, B. J.; lUya, G.; Harmandaris, V. A.; Miiller, M. M.; Kremer, K.; Deserno, M. 
Nature 2007, 447, 461-464. 

37. Blood, P. D.; Voth, G. A. Proc. Natl. Acad. Sci. USA 2006, 103, 15068-15072. 

38. Ayton, G. S.; Blood, P. D.; Voth, G. A. Biophys. J. 2007, 92, 3595-3602. 

39. Hclfrich, W. Z. Naturforsch. C 1973, 28, 693-703. 

40. Safran, S. Statistical Thermodynamics of Surfaces, Interfaces, and Membranes; Addison- 
Wesley Publishing Company: Cambridge, 1994. 

41. Risken, H. The Fokker-Planck Equation; Springer Verlag: Heidelberg, 2nd ed.; 1989. 

42. Hanggi, P.; Thomas, H. Phys. Rep. 1982, 88, 207-319. 

43. Frigo, M.; Johnson, S. G. "FFTW - for version 3.1.2", www.fftw.org, 2006. 



16 



